Emphysema represents distortion in the lung parenchyma that occurs in patients with chronic obstructive pulmonary disease (COPD) and is pathologically defined as permanent enlargement of the alveoli along with destruction of the alveolar walls. Emphysema can lead to lung hyperinflation and air trapping, which contributes to refractory dyspnea. In carefully selected patients with emphysema, lung volume reduction could improve breathlessness, lung function, quality of life and exercise capacity by decreasing lung hyperinflation 1.
Two mechanisms exist for lung volume reduction, the first, lung volume reduction surgery (LVRS), entails wedge resection of the emphysematous tissue 2. Surgical morbidity can be significant, and certain comorbidities can preclude treatment in many patients. The second is bronchoscopic lung volume reduction (BLVR) with endobronchial valves (EBV). It is a procedure where one-way valves are placed into airway segments in the most diseased lung lobe, allowing air to escape the target area but not allowing it to go back in 3. This mechanism allows the eventual deflation of the target lung lobe and the offloading diaphragmatic stretch.
Two Large clinical trials 4,5 resulted in the FDA approval of endobronchial valves for bronchoscopic lung volume reduction. Overall, respiratory complications occur in 31–35% of patients 6. The most common adverse event associated with EBV placement is pneumothorax (approximately 25% of the cases) 7,8, particularly in the immediate post operative period. Other complications include respiratory failure, acute exacerbation of COPD (AECOPD), pneumonia and hemoptysis. The frequency of these complications is significantly lower than pneumothorax in the published clinical trials.
Since FDA approval in 2018, several centers in the United States have been offering bronchoscopic lung volume reduction for patient with severe COPD. Some centers have published individual reports of complication rates 9,10. To date, there has been no study evaluating complication rates post bronchoscopic lung volume reduction on a population level. While pneumothorax is an expected complication of BLVR, acute respiratory failure can lead to dramatic instability and if progressive, can result in endobronchial valve removal.
2 Overview of the Intervention
2.1 Bronchoscopic Lung Volume Reduction (BLVR)
2.2 Surgical Lung Volume Reductuion (LVRS)
3 Objectives
To evaluate the odds of acute respiratory failure (ARF) and pneumothorax (PTX) in adult patients with a principal diagnosis of emphysema who are undergoing elective bronchoscopic lung volume reduction (BLVR), in comparison to patients with emphysema who are undergoing surgical lung volume reduction (LVRS) in the immediate post operative period
4 Research Questions
Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural acute respiratory failure in adult patients with emphysema compared to lung volume reduction surgery?
Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural pneumothorax in adult patients with emphysema compared to lung volume reduction surgery?
5 The Data
Data were obtained from the national inpatient sample (NIS); a US administrative database published under the healthcare utilization project (HCUP) by the Agency for Healthcare Research and Quality 11 through a data use agreement.
The NIS is the largest publicly available all-payer inpatient care database in the United States, containing data on more than seven million hospital stays from nearly 1,000 U.S. hospitals including non-federal public hospitals and academic medical centers. The data include discharge-level records on demographics, diagnoses, procedures, length of hospital stay (LOS), and resource utilization. The NIS uses the International Classification of Diseases, Clinical Modification 10 (ICD-10) to code for diagnoses and procedures.
The data were de-identified, not requiring an IRB approval
6 The Participants
Adult patients with a primary diagnosis of emphysema admitted electively for lung volume reduction (surgical or bronchoscopic).
The inclusion and exclusion criteria are aimed at ensuring that the patient was solely admitted for lung volume reduction, and that surgery or endobronchial valve placement was not performed for persistent air leak.
Inclusion criteria:
Adults (18 years and older)
Elective admission
Principal diagnosis of emphysema (I10-DX1)
Principal procedure for the admission is lung volume reduction, including bronchoscopic lung volume reduction (I10-PR1)
Principal procedure on the day of admission (procedure on day 0)
Exclusion criteria
Non-elective admissions
Presence of persistent air leak diagnosis (ICD10 codes)
The sample size consists of 518 adult patients with 253 (49%) patients undergoing bronchoscopic lung volume reduction and 265 patients undergoing surgical lung volume reduction (51%).
The study is a retrospective cohort study, two groups of patients with the same baseline disease (emphysema), one group is receiving the intervention of choice (BLVR), and followed during admission for the development of the outcomes after the intervention of choice.
7 The Intervention
Bronchoscopic lung volume reduction (BLVR) is the intervention of interest. 253 (49%) of the participants out of 518 in my sample have received BLVR and will be compared to surgical lung volume reduction (LVRS), in 265 (51%) of the participants.
The intervention (BLVR) will be defined as the first procedure (I10_PR1) received by the participant on day 0 of the elective admission, using ICD-10 procedure codes: 0BH38GZ, 0BH88GZ, 0BH48GZ, 0BH58GZ, 0BH68GZ, 0BH78GZ, 0BH8GZ, 0BH98GZ.
The comparison (LVRS) will be defined as the first procedure (I10_PR1) received by the participant on day 0 of the elective admission, using ICD-10 procedure codes: 0BBC0ZZ, 0BBC4ZZ, 0BBD0ZZ, 0BBD4ZZ, 0BBF0ZZ, 0BBF4ZZ, 0BBG0ZZ, 0BBG4ZZ, 0BBH0ZZ, 0BBH4ZZ, 0BBJ0ZZ, 0BBJ4ZZ, 0BBK0ZZ, 0BBK4ZZ, 0BBL0ZZ, 0BBL4ZZ, 0BBM0ZZ, 0BBM4ZZ
8 The Covariates used in the propensity score
The propensity score model used covariates that are related to patient characteristics (age, sex, race, presence of chronic respiratory failure, insurance carrier and classification of the median income for the patient’s zip code) and hospital characteristics (hospital size, location, teaching status, hospital ownership and year of the procedure)
The covariates were measured and recorded on admission and prior to the procedure.
Acute respiratory failure is a new (or increased) need for oxygen supplementation in order to maintain adequate oxygenation.
It will be defined using ICD-10 diagnosis codes J9600, J9601, J9602, J9620, J9621, J9622, J95821, J95822 and is a binary outcome (Yes/No)
9.2 Secondary outcome: Pneumothorax (PTX)
Pneumothorax
Pneumothorax is the abnormal collection of air around the lung leading to lung collapse, and possibly causing respiratory failure or in severe forms, instability and death.
This will be defined using ICD-10 diagnosis codes J930, J9311, J9312, J9383, J939, J95811, J95812 and is a binary outcome (Yes/No)
10 Importing and cleaning the data
10.1 Importing the data
The final analytical data was prepared using SAS. This included importing the national inpatient sample (NIS) into a statistical program that can handle the large size and number of original observations (roughly 21 million observations). SAS was used to properly create the variables of interest based on ICD-10 codes. A final dataset with the intended variables was then created for the final analyses.
The data was imported in the format of sas7bdat using haven package. zap_labels and clean_names were used to remove column labels and standardize names respectively into R format.
The final analytical tibble included the covariates used in building the propensity score, as well as both numeric/ factor versions of the intervention and outcomes
There were fourteen missing observations in the race variable, and nine observations in the income variable. These constituted 4.4% of the total data, and assumed to be missing at random (MAR).
Code
blvr |>tabyl(income)
income n percent valid_percent
0-25th percentile 139 0.26833977 0.2730845
26th to 50th percentile 133 0.25675676 0.2612967
51st to 75th percentile 129 0.24903475 0.2534381
76th to 100th percentile 108 0.20849421 0.2121807
<NA> 9 0.01737452 NA
Code
blvr |>tabyl(race)
race n percent valid_percent
White 421 0.81274131 0.83531746
Black 46 0.08880309 0.09126984
Other 37 0.07142857 0.07341270
<NA> 14 0.02702703 NA
Given the small sample and the desire to keep as many observations as possible, single imputation using the mice package. Fifteen imputed data sets were created and one dataset was later used as the single imputation tibble.
Code
set.seed(123)blvr_mice <-mice(blvr, m =15, printFlag =FALSE)
Code
blvr <-complete(blvr_mice, 10) |>tibble()
Checking the variables with the missing values after imputation showed preserved ratios in each level in comparison to pre-imputation
Code
blvr |>tabyl(income)
income n percent
0-25th percentile 142 0.2741313
26th to 50th percentile 136 0.2625483
51st to 75th percentile 131 0.2528958
76th to 100th percentile 109 0.2104247
Code
blvr |>tabyl(race)
race n percent
White 432 0.83397683
Black 48 0.09266409
Other 38 0.07335907
10.5 Adding an ID row
An ID row was added to identify observations if needed
Code
blvr <- tibble::rowid_to_column(blvr, "ID"); blvr <- blvr |>mutate (ID =as.character(ID))
10.6 Saving final tibbe
Code
write.csv(blvr, "blvr.csv")
11 Table (1) And Codebook
11.1 Final Codebook
Variable
Type
Description
Role
age
Numeric
Age in years
Co-variate
sex
Factor (w/2 levels)
Sex of the subject: Male or Female
Covariate
race
Factor (w/3 levels)
White, Black, Other
Covariate
crf
Factor (w/2 levels)
Presence of chronic respiratory failure: Yes or No
Covariate
income
Factor (w/4 levels)
Quartile classification of the estimated median household income of residents in the patient’s ZIP Code: 0-25th percentile, 26-50th percentile, 51-75th percentile and 76-00th percentile
Covariate
insurance
Factor (w/3 levels)
Medicare, Medicaid or Private
Covariate
hospital_size
Factor (w/3 levels)
Small, Medium, Large
Covariate
teaching
Factor (w/3 levels)
Teaching status of the hospital: Rural, Urban non-teaching and Urban teaching
Covariate
region
Factor (w/4 levels)
Region of the hospital: Northeast, Midwest, South and West
Covariate
ownership
Factor (w/3 levels)
Ownership of the hospital: Government, Private, non-profit and Private for-profit
Covariate
year
Character
Year of the procedure
Covariate
vol_reduction
numeric
Type of lung volume reduction: bronchoscopic (1) or surgical (0)
Intervention/Comparator
vol_reduction_f
Factor (w/2 levels)
Type of lung volume reduction: bronchoscopic BLVR or surgical LVRS
Based on table 1, patients receiving bronchoscopic lung volume reduction (BLVR) were noted to be older (mean age of 67 years), with a higher percentage of females (43.5%), likely due to a preference of a less invasive approach compared to the surgical group. BLVR patients also tended to have a higher percentage of chronic respiratory failure, likely due to the ease of valve removal should the patient’s baseline respiratory failure worsen after procedure. The patient population was more white (90%), with medicare as the primary insurance provider (77.5%) which is likely due to the age bracket in this group (65 years and older).
As expected both groups had a similar distribution of median income per zip code, and majority of these procedures were done in large, non for profit and urban teaching hospitals. Lung volume reduction surgery number distribution was steady over the years included (roughly around 30%), however bronchoscopic lung volume reduction numbers significantly increased over the years (starting in 2019 and 2020)
12 Unadjusted Outcomes
12.1 Acute Respiratory Failure arf
Code
ggplot(blvr, aes(x = vol_reduction_f, fill = arf_f)) +geom_bar(aes(y = (after_stat(count))/sum(after_stat(count))), position ="dodge") +scale_y_continuous(labels = scales::percent) +scale_fill_viridis_d(option ="turbo", alpha =0.5) +labs(title ="Respiratory Failure", x ="Volume Reduction", y ="% of Patients", fill ="Resp. Failure")
The visualization above suggested that BLVR patients tended to suffer a higher percentage of acute respiratory failure
A logistic regression model was fit to assess the odds of respiratory failure based on the type of volume reduction received, then the coefficient was exponentiated to obtain an odds ratio
If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 1.65 odds of acute respiratory failure (95% CI (1,2.74)). This confidence interval contained 1, thus it was difficult to reliably say that BLVR would meaningfully lead to a higher odds of acute respiratory failure
Overall, it seems that the patients receiving BLVR tended to have higher odds of respiratory failure, even though the odds ratio’s interval contained 1, this odds ratio shows an interesting trend
12.2 Pneumothorax ptx
Code
ggplot(blvr, aes(x = vol_reduction_f, fill = ptx_f)) +geom_bar(aes(y = (after_stat(count))/sum(after_stat(count))), position ="dodge") +scale_y_continuous(labels = scales::percent) +scale_fill_viridis_d(option ="turbo", alpha =0.5) +labs(title ="Pneumothorax", x ="Volume Reduction", y ="% of Patients", fill ="PTX")
The visualization above suggested that BLVR patients tended to suffer a lower percentage of pneumothorax
A logistic regression model was fit to assess the odds of pneumothorax based on the type of volume reduction received, then the coefficient was exponentiated to obtain an odds ratio
If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR had a 0.4 odds of pneumothorax (95% CI (0.27,0.58)), which is lower odds than the subject receiving LVRS, which is plausible as LVRS is more invasive.
13 Propensity Score
13.1 Fitting the model
A logistic regression model was used to estimate the propensity score. Having BLVR/LVRS was modeled as the outcome, with patient and hospital characteristics as independent variables. Patient characteristics used in the propensity score included age, sex, race, crf, income and insurance. Hospital characteristics included hospital_size, teaching status, region, ownership and year of the procedure. The raw and linear propensity score values were saved to the blvr tibble.
Code
psmodel <-glm(vol_reduction ~ age + sex + race + crf + income + insurance + hospital_size + teaching + region + ownership + year,family=binomial, data=blvr)summary(psmodel)
Call:
glm(formula = vol_reduction ~ age + sex + race + crf + income +
insurance + hospital_size + teaching + region + ownership +
year, family = binomial, data = blvr)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -4.696344 1.503545 -3.124 0.001787 **
age 0.064113 0.015710 4.081 4.48e-05 ***
sexFemale 0.716995 0.258418 2.775 0.005528 **
raceBlack -0.481109 0.451914 -1.065 0.287055
raceOther -1.521225 0.599723 -2.537 0.011195 *
crfYes 1.294080 0.288002 4.493 7.01e-06 ***
income26th to 50th percentile 0.393352 0.350994 1.121 0.262423
income51st to 75th percentile -0.286941 0.342652 -0.837 0.402361
income76th to 100th percentile 0.026792 0.386403 0.069 0.944722
insuranceMedicaid -1.353403 0.515379 -2.626 0.008639 **
insurancePrivate -0.469448 0.318487 -1.474 0.140484
hospital_sizeMedium -2.400828 0.671416 -3.576 0.000349 ***
hospital_sizeLarge -2.077553 0.616456 -3.370 0.000751 ***
teachingUrban nonteaching 0.678464 0.954277 0.711 0.477102
teachingUrban teaching -0.472681 0.823311 -0.574 0.565885
regionMidwest 0.044408 0.388618 0.114 0.909021
regionSouth 0.668451 0.378964 1.764 0.077750 .
regionWest 0.167103 0.437435 0.382 0.702457
ownershipPrivate, non-profit -0.006272 0.357697 -0.018 0.986011
ownershipPrivate, for-profit -1.513048 0.774833 -1.953 0.050850 .
year2019 2.211839 0.517582 4.273 1.93e-05 ***
year2020 3.281191 0.507832 6.461 1.04e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 717.82 on 517 degrees of freedom
Residual deviance: 436.37 on 496 degrees of freedom
AIC: 480.37
Number of Fisher Scoring iterations: 6
ggplot(blvr, aes(x = vol_reduction_f, y = ps, color = vol_reduction_f)) +geom_boxplot() +geom_jitter(width =0.1, alpha =0.3) +scale_color_viridis_d(option ="plasma", end =0.5) +labs(title ="Boxplot & Jitter", subtitle ="Distribution Of the Raw PS", x ="Type of Procedure", y ="Propensity Scores",caption ="* Line indicates sample median", col ="Procedure")
Code
ggplot(blvr, aes(x = linps, fill = vol_reduction_f)) +geom_density(alpha =0.3) +scale_fill_viridis_d(option ="turbo") +labs(title ="Density Plot", subtitle ="For the Linear Propensity Score", x ="Linps", y ="Density",fill ="Volume Reduction")
The overlap between the propensity scores in intervention groups was assessed numerically and graphically. The boxplot and the density plots show some overlap in the propensity scores between the intervention and the treatment.
The LVRS group had a lower propensity score mean, with problematic values falling below 0.01. The BLVR group also had some problematic high values over 0.99. Troubleshooting will be addressed at a later step
Rubin’s rule 1 value was higher than 50%, this meant that running an unadjusted model could not be justified without applying the propensity scores to account for selection bias.
Rubin’s rule 2 was on the lower bound of acceptable value as well, as the variance ratio of the linear propensity score in the BLVR (intervention) group to the variance of the linear propensity score in the LVRS (control) group should be close to 1 (and ideally between 0.8 and 1.25)
13.4 Troubleshooting the Extreme Propensity Score Values
Code
blvr |>count(ps >=0.99); blvr |>count(ps <0.01)
# A tibble: 2 × 2
`ps >= 0.99` n
<lgl> <int>
1 FALSE 514
2 TRUE 4
# A tibble: 2 × 2
`ps < 0.01` n
<lgl> <int>
1 FALSE 490
2 TRUE 28
There were twenty eight observations with propensity score values below 0.01 and four observations with propensity score values above 0.99. The extreme lows seemed to be present in the LVRS group while the extreme highs were present in the BLVR group
Code
mosaic::favstats(ps ~ vol_reduction_f, data = blvr)
Registered S3 method overwritten by 'mosaic':
method from
fortify.SpatialPolygonsDataFrame ggplot2
vol_reduction_f min Q1 median Q3 max
1 LVRS 0.0005808267 0.03308973 0.1778799 0.4680281 0.9256507
2 BLVR 0.0263071482 0.59915004 0.7983214 0.9061054 0.9964521
mean sd n missing
1 0.2679640 0.2584781 265 0
2 0.7193262 0.2402497 253 0
The following strategies were attempted in order to resolve the above issue:
Colinearity check using vif didn’t reveal problematic colinearity between the variables used in the propensity score model
Code
car::vif(psmodel)
GVIF Df GVIF^(1/(2*Df))
age 1.634643 1 1.278532
sex 1.070972 1 1.034878
race 1.204488 2 1.047612
crf 1.105123 1 1.051248
income 1.428943 3 1.061294
insurance 1.667039 2 1.136283
hospital_size 1.403948 2 1.088523
teaching 1.250028 2 1.057377
region 1.470981 3 1.066435
ownership 1.355067 2 1.078922
year 1.164515 2 1.038810
Re-constructing the propensity score model, adding one covariate at a time, while checking which variable drives the propensity score values below 0.1 or above 0.99 showed that each of the following covariates race, year, sex, ownership and crf could independently drop the propensity score below 0.01 or drive it up above 0.99. Excluding these covariates form the model would be problematic in itself
Collapsing multicategorical variables was attempted, these variables included year, race, ownership, insurance, hospital_size, teaching and region. It did not result in the improvement of the propensity score even when each covariate was added individually after collapsing.
Filtering obaservations with problematic propensity scores proved to be the most reasonable approach in this situation. This resulted in a sample with 237 observations on the control intervention (LVRS) and 249 observation on the studied intervention (BLVR). The resulting dataset blvr_matching was used for propensity score matching analyses
mosaic::favstats(ps ~ vol_reduction_f, data = blvr_matching)
vol_reduction_f min Q1 median Q3 max mean
1 LVRS 0.01020130 0.05651206 0.2244966 0.4944178 0.9256507 0.2991748
2 BLVR 0.02630715 0.58669227 0.7927773 0.9014236 0.9892394 0.7149098
sd n missing
1 0.2558878 237 0
2 0.2396084 249 0
A density plot was created to describe the overlap of the linear propensity scores between the two intervention groups after balancing the dataset, the overlap is described graphically in the density plot below.
Code
ggplot(blvr_matching, aes(x = linps, fill = vol_reduction_f)) +geom_density(alpha =0.3) +scale_fill_viridis_d(option ="turbo") +labs(title ="Density Plot", subtitle ="For the Linear Propensity Score", x ="Linps", y ="Density",fill ="Volume Reduction")
14 Approach (1): Matching
Propensity score matching was conducted using a 1:1 caliper matching on the linear propensity score without replacement from the Matching package. A default caliper value of 0.2 was selected. The type of matching was chosen based on the similar sample size between the intervention and control groups.
Code
X <- blvr_matching$linpsTr <-as.logical(blvr_matching$vol_reduction)match <-Match(Tr = Tr, X = X, M =1, caliper =0.2,estimand ="ATT", replace =FALSE, ties =FALSE)
Warning in Match(Tr = Tr, X = X, M = 1, caliper = 0.2, estimand = "ATT", :
replace==FALSE, but there are more (weighted) treated obs than control obs.
Some treated obs will not be matched. You may want to estimate ATC instead.
Code
summary(match)
Estimate... 0
SE......... 0
T-stat..... NaN
p.val...... NA
Original number of observations.............. 486
Original number of treated obs............... 249
Matched number of observations............... 107
Matched number of observations (unweighted). 107
Caliper (SDs)........................................ 0.2
Number of obs dropped by 'exact' or 'caliper' 142
One hundred and seven BLVR patients were matched to one hundred and seven LVRS patients. The degree of covariate imbalance was assessed before and after propensity score matching for each covariate as well as the linear and raw propensity scores
bal.plot(match,treat = blvr_matching$vol_reduction,covs = covs,var.name ="linps",which ="both",sample.names =c("Unmatched Sample", "Matched Sample")) +scale_fill_viridis_d(option ="turbo") +labs(title ="Distributional Balance for linear PS", subtitle ="1:1 Caliper Matching", x ="Linear PS")
A Love plot of standardized differences before and after matching was obtained to evaluate matching success, showing improvement over the unmatched sample
If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 1.73 odds of acute respiratory failure (95% CI (0.82,3.63)). The subject receiving BLVR tended to have a higher odds ratio of respiratory failue in comparison to patients receiving LVRS, however this confidence interval crosses 1
If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 0.68 odds of pneumothorax (95% CI (0.38,1.22)). Patients receiving BLVR seemed to have a higher odds ratio of pneumothorax in comparison to patients receiving BLVR, however this confidence interval contained 1, the results were not statistically detectable and a stability analysis instead of a sensitivity analysis was carried out.
15.2.1 Stability Analysis
Trying different imputed datasets (generated through mice package) resulted in similar outcome estimates (odds ratio)
The matching strategies described below were attempted. None of these strategies provided an improvement in the effective sample size, propensity score distributional balance, standardized differences or Rubin’s rules in comparison to the matching strategy used in this analysis (1:1 matching with a caliper)
Match
Trtmnt
Ctrl
Distributional Balance
Std Differences
Rubin’s 1
Rubin’s 2
1:1 Greedy W replacement
249
75
Better overlap, not ideal
better but dispersed
8.8
1.27
Genetic matching
273
78
Better overlap, not ideal
Dispersed
8.8
1.26
1:1 Neartest neighbor matching W Replacement
249
79
Better overlap, not ideal
better but dispersed (> 0.1)
8
1.21
The following steps attempted comparing the odds ratios obtained through matching to the odds ratios after ATT weighting and a double robust approach (by adjusting for the linear propensity score in the regression model)
16 Approach (2): Propensity Score Weighting
Given the initial difference in sample sizes across the two volume reduction cohorts, propensity score matching has left a selection of the patients unmatched. Weighting by the inverse propensity score (average of the treatment effect on the treated) was next carried out using the initial database blvr (the database that still contained propensity scores below 0.01 and above 0.99)
Then the resulting ATT (average treatment effect on the treated) weights were plotted in the figure below, where the treatment group got the same weight of 1, and the control group had varying weights
Code
ggplot(blvr, aes(x = ps, y = wts1, color = vol_reduction_f)) +geom_point(size =2, alpha =0.4) +facet_wrap(~ vol_reduction_f) +labs(x ="Estimated Propensity for Volume Reduction",y ="ATT weights for blvr",title ="ATT weighting structure",subtitle ="In the BLVR dataset",col ="Volume Reduction")
The balance after weighting was assessed using dx.wts from the twang package. Prior to weighting, the control sample was 237 participants, and after weighting, the effective sample size was 54 participants, representing a sample size that would be required to achieve the same level of precision if that sample were a simple random sample
type n.treat n.ctrl ess.treat ess.ctrl max.es mean.es max.ks
1 unw 253 265 253 265.0000 2.0306797 0.4589984 0.6194496
2 253 265 253 54.7126 0.3077257 0.1118164 0.2134387
mean.ks iter
1 0.14873095 NA
2 0.05053566 NA
The standardized differences were compared before and after weighting using a Love.plot. The weighting process showed a more satisfactory balance of standardized differences in comparison to the unweighted standardized differences.
Code
bal.table(bal.wts1)
$unw
tx.mn tx.sd ct.mn ct.sd std.eff.sz stat
age 67.881 8.390 56.400 15.065 1.368 10.789
sex:Male 0.565 0.496 0.713 0.452 -0.299 12.293
sex:Female 0.435 0.496 0.287 0.452 0.299 NA
race:White 0.901 0.298 0.770 0.421 0.440 8.925
race:Black 0.067 0.250 0.117 0.321 -0.199 NA
race:Other 0.032 0.175 0.113 0.317 -0.466 NA
crf:No 0.538 0.499 0.872 0.334 -0.670 69.838
crf:Yes 0.462 0.499 0.128 0.334 0.670 NA
insurance:Medicare 0.775 0.418 0.408 0.491 0.879 37.291
insurance:Medicaid 0.043 0.204 0.192 0.394 -0.731 NA
insurance:Private 0.182 0.386 0.400 0.490 -0.566 NA
income:0-25th percentile 0.269 0.443 0.279 0.449 -0.024 0.225
income:26th to 50th percentile 0.269 0.443 0.257 0.437 0.027 NA
income:51st to 75th percentile 0.241 0.428 0.264 0.441 -0.054 NA
income:76th to 100th percentile 0.221 0.415 0.200 0.400 0.051 NA
hospital_size:Small 0.107 0.309 0.045 0.208 0.199 3.868
hospital_size:Medium 0.166 0.372 0.208 0.406 -0.112 NA
hospital_size:Large 0.727 0.445 0.747 0.435 -0.045 NA
teaching:Rural 0.028 0.164 0.026 0.160 0.008 0.965
teaching:Urban nonteaching 0.111 0.314 0.075 0.264 0.112 NA
teaching:Urban teaching 0.862 0.345 0.898 0.302 -0.106 NA
region:Northeast 0.126 0.332 0.238 0.426 -0.335 4.752
region:Midwest 0.304 0.460 0.298 0.457 0.014 NA
region:South 0.415 0.493 0.294 0.456 0.245 NA
region:West 0.154 0.361 0.170 0.375 -0.043 NA
ownership:Government 0.123 0.328 0.140 0.347 -0.052 2.065
ownership:Private, non-profit 0.846 0.361 0.792 0.406 0.148 NA
ownership:Private, for-profit 0.032 0.175 0.068 0.252 -0.207 NA
year:2018 0.024 0.152 0.321 0.467 -1.952 50.600
year:2019 0.300 0.458 0.366 0.482 -0.143 NA
year:2020 0.676 0.468 0.313 0.464 0.775 NA
ps 0.719 0.240 0.268 0.258 1.879 20.616
linps 1.288 1.526 -1.811 2.145 2.031 19.031
p ks ks.pval
age 0.000 0.359 0.000
sex:Male 0.000 0.148 0.000
sex:Female NA 0.148 0.000
race:White 0.000 0.131 0.000
race:Black NA 0.050 0.000
race:Other NA 0.082 0.000
crf:No 0.000 0.334 0.000
crf:Yes NA 0.334 0.000
insurance:Medicare 0.000 0.367 0.000
insurance:Medicaid NA 0.149 0.000
insurance:Private NA 0.218 0.000
income:0-25th percentile 0.879 0.010 0.879
income:26th to 50th percentile NA 0.012 0.879
income:51st to 75th percentile NA 0.023 0.879
income:76th to 100th percentile NA 0.021 0.879
hospital_size:Small 0.021 0.061 0.021
hospital_size:Medium NA 0.042 0.021
hospital_size:Large NA 0.020 0.021
teaching:Rural 0.381 0.001 0.381
teaching:Urban nonteaching NA 0.035 0.381
teaching:Urban teaching NA 0.036 0.381
region:Northeast 0.003 0.111 0.003
region:Midwest NA 0.006 0.003
region:South NA 0.121 0.003
region:West NA 0.016 0.003
ownership:Government 0.127 0.017 0.127
ownership:Private, non-profit NA 0.053 0.127
ownership:Private, for-profit NA 0.036 0.127
year:2018 0.000 0.297 0.000
year:2019 NA 0.066 0.000
year:2020 NA 0.363 0.000
ps 0.000 0.619 0.000
linps 0.000 0.619 0.000
[[2]]
tx.mn tx.sd ct.mn ct.sd std.eff.sz stat
age 67.881 8.390 66.767 8.471 0.133 0.923
sex:Male 0.565 0.496 0.613 0.487 -0.096 0.372
sex:Female 0.435 0.496 0.387 0.487 0.096 NA
race:White 0.901 0.298 0.866 0.341 0.118 0.576
race:Black 0.067 0.250 0.104 0.306 -0.148 NA
race:Other 0.032 0.175 0.030 0.170 0.010 NA
crf:No 0.538 0.499 0.655 0.475 -0.236 2.117
crf:Yes 0.462 0.499 0.345 0.475 0.236 NA
insurance:Medicare 0.775 0.418 0.782 0.413 -0.018 0.111
insurance:Medicaid 0.043 0.204 0.034 0.182 0.046 NA
insurance:Private 0.182 0.386 0.184 0.387 -0.005 NA
income:0-25th percentile 0.269 0.443 0.325 0.469 -0.128 0.790
income:26th to 50th percentile 0.269 0.443 0.180 0.384 0.200 NA
income:51st to 75th percentile 0.241 0.428 0.265 0.442 -0.057 NA
income:76th to 100th percentile 0.221 0.415 0.229 0.420 -0.019 NA
hospital_size:Small 0.107 0.309 0.101 0.301 0.019 0.006
hospital_size:Medium 0.166 0.372 0.169 0.375 -0.009 NA
hospital_size:Large 0.727 0.445 0.730 0.444 -0.006 NA
teaching:Rural 0.028 0.164 0.017 0.130 0.063 3.263
teaching:Urban nonteaching 0.111 0.314 0.037 0.188 0.236 NA
teaching:Urban teaching 0.862 0.345 0.946 0.226 -0.245 NA
region:Northeast 0.126 0.332 0.208 0.406 -0.247 1.450
region:Midwest 0.304 0.460 0.334 0.472 -0.065 NA
region:South 0.415 0.493 0.363 0.481 0.105 NA
region:West 0.154 0.361 0.094 0.292 0.166 NA
ownership:Government 0.123 0.328 0.099 0.299 0.071 1.134
ownership:Private, non-profit 0.846 0.361 0.887 0.316 -0.115 NA
ownership:Private, for-profit 0.032 0.175 0.013 0.115 0.104 NA
year:2018 0.024 0.152 0.032 0.176 -0.055 0.187
year:2019 0.300 0.458 0.318 0.466 -0.039 NA
year:2020 0.676 0.468 0.650 0.477 0.056 NA
ps 0.719 0.240 0.663 0.223 0.233 1.767
linps 1.288 1.526 0.818 1.184 0.308 2.489
p ks ks.pval
age 0.356 0.099 0.734
sex:Male 0.542 0.048 0.542
sex:Female NA 0.048 0.542
race:White 0.552 0.035 0.552
race:Black NA 0.037 0.552
race:Other NA 0.002 0.552
crf:No 0.146 0.118 0.146
crf:Yes NA 0.118 0.146
insurance:Medicare 0.870 0.008 0.870
insurance:Medicaid NA 0.009 0.870
insurance:Private NA 0.002 0.870
income:0-25th percentile 0.494 0.057 0.494
income:26th to 50th percentile NA 0.089 0.494
income:51st to 75th percentile NA 0.024 0.494
income:76th to 100th percentile NA 0.008 0.494
hospital_size:Small 0.990 0.006 0.990
hospital_size:Medium NA 0.003 0.990
hospital_size:Large NA 0.003 0.990
teaching:Rural 0.044 0.010 0.044
teaching:Urban nonteaching NA 0.074 0.044
teaching:Urban teaching NA 0.084 0.044
region:Northeast 0.229 0.082 0.229
region:Midwest NA 0.030 0.229
region:South NA 0.052 0.229
region:West NA 0.060 0.229
ownership:Government 0.319 0.023 0.319
ownership:Private, non-profit NA 0.042 0.319
ownership:Private, for-profit NA 0.018 0.319
year:2018 0.784 0.008 0.784
year:2019 NA 0.018 0.784
year:2020 NA 0.026 0.784
ps 0.078 0.213 0.029
linps 0.013 0.213 0.029
ggplot(balance.att.weights, aes(x = szd, y =reorder(names, szd), color = timing)) +geom_point(size =2, alpha =0.5) +geom_vline(xintercept =0) +geom_vline(xintercept =c(-10,10), linetype ="dashed", col ="blue") +labs(x ="Standardized Differences", y ="", title ="Love Plot of Standarized Diffs in The BLVR Data",subtitle ="Before and After ATT Weighting")
Rubin’s rules were obtained after weighting. Rubin’s rule 1 was obtained by using the standardized effect size after weighting for the linear propensity score (0.28). When multiplied by 100, the result was 28% indicating that Rubin’s rule 1 was satisfied (although not perfectly). Rubin’s rule 2 was calculated by dividing the square root of treatment group standard deviation by the square root of control group standard deviation, resulting in a value of 1.54. Rubin’s rule 2 was also satisfied, although less perfectly.
The estimate of treatment effect on acute respiratory failure was obtained by a regression model using survey weighting and a quasi binomial design.
The weighted odds ratio estimate for the BLVR group was 1.68, with a 95% confidence interval (0.83,3.4). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having acute respiratory failure for the person receiving BLVR was higher and 1.68 the odds of respiratory failure compared to the person receiving LVRS. It should be noted that this confidence interveal crosses 1
The weighted odds ratio estimate for the `BLVR group was 0.36, with a 95% confidence interval (0.19,0.67). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having pneumothorax for the person receiving BLVR was lower and was 0.36 the odds of pneumothorax of the person receiving LVRS. The confidence interval did not cross 1
18 Apprach (3): A Double-Robust Estimate
A logistic regression for the outcomes of acute respiratory failure and pneumothorax was fitted with a regression adjustment for the linear propensity score to obtain a doubly-robust estimate and account for residual imbalance after weighting
The doubly-robust odds ratio estimate for the BLVR group was 1.89, with a 95% confidence interval (0.93,3.87). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having acute respiratory failure for the person receiving BLVR was higher and 1.89 the odds of acute respiratory failure compared to the person receiving LVRS. The confidence interval here crosses 1
The doubly-robust odds ratio estimate for the `BLVR group was 0.38, with a 95% confidence interval (0.2,0.75). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having pneumothorax for the person receiving BLVR was lower and was 0.38 the odds of pneumothorax of the person receiving LVRS. The confidence interval did not cross 1
19 Final comparison of Outcome Estimates
19.1 Acute Respiratory Failure arf
Code
arf_blvr <-tibble(analysis =c("Unadjusted", "Matched", "ATT Weighted", "Double-Robust"),estimate =c(1.65, 1.73, 1.68, 1.89),conf.low =c(1, 0.82, 0.83, 0.93),conf.high =c(2.74, 3.63, 3.4, 3.87))ggplot(arf_blvr, aes(x = analysis, y = estimate)) +geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width =0.5, color ="orange3") +geom_label(aes(label = estimate), size =6, color ="white", fill ="lightblue3") +labs(title ="Comparison Of Odds Ratio Estimates", subtitle ="Between Different Adjustment Approaches In Respiratory Failure", x ="Method of PS Adjustment", y ="Odds Ratio")
Method
Estimate
Lower 95% CI
upper 95% CI
Unadjusted
1.65
1.00
2.74
Matched
1.73
0.82
3.63
Weighted
1.68
0.83
3.4
Double-Robust
1.89
0.93
3.87
Based on the comparison above, the odds ratio of acute respiratory failure didn’t substantially change from its estimate prior to propensity score adjustment. All of the confidence intervals contain 1.
19.2 Pneumothorax ptx
Code
ptx_blvr <-tibble(analysis =c("Unadjusted", "Matched", "ATT Weighted", "Double-Robust"),estimate =c(0.4, 0.68, 0.36, 0.38),conf.low =c(0.27, 0.38, 0.19, 0.2),conf.high =c(0.58, 1.22, 0.67, 0.75))ggplot(ptx_blvr, aes(x = analysis, y = estimate)) +geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width =0.5, color ="purple3") +geom_label(aes(label = estimate), size =6, color ="white", fill ="lightpink3") +labs(title ="Comparison Of Odds Ratio Estimates", subtitle ="Between Different Adjustment Approaches In Pneumothorax", x ="Method of PS Adjustment", y ="Odds Ratio")
Method
Estimate
Lower 95% CI
upper 95% CI
Unadjusted
0.4
0.27
0.58
Matched
0.68
0.38
1.22
Weighted
0.36
0.19
0.67
Double-Robust
0.38
0.2
0.75
The comparison above showed that the odds ratio for pneumothorax also didn’t substantially change after propensity score adjustment, reflecting a plausible conclusion that patients undergoing lung volume reduction surgery tended to have higher odds of pneumothorax based on the nature of invasive surgical intervention on the lung.
20 Consclusions
20.1 Answering The Research Questions
Is bronchoscopic lung volume reduction associated with higher odds of developing post-procedural acute respiratory failure in adult patients with emphysema compared to lung volume reduction surgery?
Yes. Based on the odds ratio obtained prior to propensity score adjustment and after, the patients undergoing bronchoscopic lung volume reduction surgery tended to have higher odds of developing acute respiratory failure, granted that the confidence interval for this estimate always crossed 1
It was noted that despite propensity score matching, weighting and a double robust adjustment (Section 19.1), the odds ratio of acute respiratory failure did not show a meaningful change from prior to adjustment. The small fluctuations wouldn’t have significant clinical implication.
Even though the confidence intervals crossed 1, this analysis showed an interesting trend in the odds of respiratory failure. A possible explanation could be that there was a higher percentage of patients with chronic respiratory failure undergoing bronchoscopic lung volume reduction in comparison to patients undergoing surgical lung volume reduction, so the patient’s lungs have already been compromised and further interventions
Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural pneumothorax in adult patients with emphysema compared to lung volume reduction surgery?
No. BLVR was associated with lower odds of pneumothorax (Section 19.2) despite propensity score adjustment methods, which was rather plausible, as the surgical procedure was more invasive and expected to cause a higher incidence of pneumothorax. The presence of pneumothorax is well described in the literature, and could indicate a favorable response after lung volume reduction
20.2 Future Directions
Code
ggplot(blvr, aes(x = year, fill = vol_reduction_f)) +geom_bar(position ="dodge", col ="white") +scale_fill_viridis_d(option ="turbo", alpha =0.5) +labs(title ="Lung Volume Reduction", subtitle ="Distribution By Year", x ="Volume Reduction", y ="# of Patients", fill ="Volume Reduction")
Bronchoscopic lung volume reduction with endobronchial valves has largely replaced lung volume reduction surgery as a possible solution to improve excessive dyspnea. Patient selection is a key element in developing post procedural complications that could be life threatening.
Understanding the potential complications that may arise in real-world settings outside controlled clinical trials becomes increasingly crucial as the number of procedures grows and their availability expands beyond large referral centers. This knowledge aids physicians in better patient selection for the procedure, enabling them to implement effective action plans. For instance, the introduction of a “Pneumothorax Travel Cart” containing procedural kits for pneumothorax management has significantly improved time to management of pneumothorax after bronchoscopic lung volume reduction.
By trying to get a sense of the complication rates nationally, this information could be used to obtain further data from endobronchial valve device companies. These data could include more patient covariates that play a role in procedural selection, and provide a better characterization of the complications encountered.
20.3 Learning Curve
By doing this project, I got a chance to apply all the matching methods (Section 14) we learned over the span of this course. While I was slightly distraught at the beginning when I ended up with a sample with more treated patients than controls, I used it as an opportunity to go over all the plausible methods to properly troubleshoot and be prepared for similar scenarios in the future beyond the security of this class.
21 Acknowledgments
Dr. Marc Georgi, a gastroenterologist at Inova Fairfax Hospital, for simplifying the use of the Healthcare Utilization Project databases (the national inpatient sample and the national readmission database).
Dr. Thomas Love, for making statistics incredibly enjoyable!
22 References
Global strategy for the diagnosis, management and prevention of chronic obstructive pulmonary disease. Global Initiative for Chronic Obstructive Lung Disease; 2021
Fishman A, Martinez F, Naunheim K, et al. A randomized trial comparing lung-volume-reduction surgery with medical therapy for severe emphysema. N Engl J Med. 2003;348(21):2059-2073
Fernandez-Bussy S, Labarca G, Herth FJF. Bronchoscopic Lung Volume Reduction in Patients with Severe Emphysema. Semin Respir Crit Care Med. 2018;39(6):685-692.
Criner GJ, Sue R, Wright S, et al. A Multicenter Randomized Controlled Trial of Zephyr Endobronchial Valve Treatment in Heterogeneous Emphysema (LIBERATE). Am J Respir Crit Care Med. 2018;198(9):1151-1164.
Criner GJ, Delage A, Voelker K, et al. Improving Lung Function in Severe Heterogenous Emphysema with the Spiration Valve System (EMPROVE). A Multicenter, Open-Label Randomized Controlled Clinical Trial. Am J Respir Crit Care Med. 2019;200(11):1354-1362
Zantah M, Gangemi AJ, Criner GJ. Bronchoscopic lung volume reduction: status quo. Annals of translational medicine 2020; 8:1469.
Davey C, Zoumot Z, Jordan S, McNulty WH, Carr DH, Hind MD, Hansell DM, Rubens MB, Banya W, Polkey MI, et al. Bronchoscopic lung volume reduction with endobronchial valves for patients with heterogeneous emphysema and intact interlobar fissures (the BeLieVeR-HIFi study): a randomised controlled trial. Lancet 2015;386:1066–1073
Valipour A, Slebos DJ, Herth F, Darwiche K, Wagner M, Ficker JH, Petermann C, Hubner RH, Stanzel F, Eberhardt R; IMPACT Study Team. Endobronchial valve therapy in patients with homogeneous emphysema: results from the IMPACT study. Am J Respir Crit Care Med 2016;194:1073–1082.
Fiorelli A, D’Andrilli A, Bezzi M, et al. Complications related to endoscopic lung volume reduction for emphysema with endobronchial valves: results of a multicenter study. J Thorac Dis. 2018;10(Suppl 27):S3315-S3325. doi:10.21037/jtd.2018.06.69
Posthuma R, Vaes AW, Walraven KHM, et al. Implementation of Bronchoscopic Lung Volume Reduction Using One-Way Endobronchial Valves: A Retrospective Single-Centre Cohort Study. Respiration. 2022;101(5):476-484. doi:10.1159/000520885
Healthcare Cost and Utilization Project. HCUP National Inpatient Sample (NIS). Agency for Healthcare Research and Quality; 2019-2020.
---title: "Complications in Bronchoscopic Lung Volume Reduction"subtitle: "Insights form the National Inpatient Sample"author: "Hala Nas"date: last-modifiedformat: html: toc: true number-sections: true code-fold: show code-tools: true code-overflow: wrap embed-resources: true date-format: iso theme: spacelab---# R Setup {.unnumbered}```{r}#| message: false#| warning: falseknitr::opts_chunk$set(comment =NA) library(janitor)library(broom)library(haven)library(naniar)library(mice)library(kableExtra)library(tableone)library(skimr)library(Matching)library(cobalt)library(survival)library(twang)library(survey)library(conflicted)library(knitr)library(tidyverse)conflicts_prefer(dplyr::select(), dplyr::filter())theme_set(theme_bw()) ```# Background Emphysema represents distortion in the lung parenchyma that occurs in patients with chronic obstructive pulmonary disease (COPD) and is pathologically defined as permanent enlargement of the alveoli along with destruction of the alveolar walls. Emphysema can lead to lung hyperinflation and air trapping, which contributes to refractory dyspnea. In carefully selected patients with emphysema, lung volume reduction could improve breathlessness, lung function, quality of life and exercise capacity by decreasing lung hyperinflation ^1^. Two mechanisms exist for lung volume reduction, the first, lung volume reduction surgery (LVRS), entails wedge resection of the emphysematous tissue ^2^. Surgical morbidity can be significant, and certain comorbidities can preclude treatment in many patients. The second is bronchoscopic lung volume reduction (BLVR) with endobronchial valves (EBV). It is a procedure where one-way valves are placed into airway segments in the most diseased lung lobe, allowing air to escape the target area but not allowing it to go back in ^3^. This mechanism allows the eventual deflation of the target lung lobe and the offloading diaphragmatic stretch.Two Large clinical trials ^4^^,^^5^ resulted in the FDA approval of endobronchial valves for bronchoscopic lung volume reduction. Overall, respiratory complications occur in 31–35% of patients ^6^. The most common adverse event associated with EBV placement is pneumothorax (approximately 25% of the cases) ^7^^,^^8^, particularly in the immediate post operative period. Other complications include respiratory failure, acute exacerbation of COPD (AECOPD), pneumonia and hemoptysis. The frequency of these complications is significantly lower than pneumothorax in the published clinical trials.Since FDA approval in 2018, several centers in the United States have been offering bronchoscopic lung volume reduction for patient with severe COPD. Some centers have published individual reports of complication rates ^9^^,^^10^. To date, there has been no study evaluating complication rates post bronchoscopic lung volume reduction on a population level. While pneumothorax is an expected complication of BLVR, acute respiratory failure can lead to dramatic instability and if progressive, can result in endobronchial valve removal. # Overview of the Intervention ## Bronchoscopic Lung Volume Reduction (BLVR){{< video https://youtu.be/KFHTxXLQqug?feature=shared >}}## Surgical Lung Volume Reductuion (LVRS){{< video https://youtu.be/MJL_wEuk6zc?feature=shared >}}# ObjectivesTo evaluate the odds of acute respiratory failure (ARF) and pneumothorax (PTX) in adult patients with a principal diagnosis of **emphysema** who are undergoing **elective** bronchoscopic lung volume reduction (BLVR), in comparison to patients with **emphysema** who are undergoing surgical lung volume reduction (LVRS) in the immediate post operative period # Research Questions> Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural **acute respiratory failure** in adult patients with emphysema compared to lung volume reduction surgery?> Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural **pneumothorax** in adult patients with emphysema compared to lung volume reduction surgery?# The Data Data were obtained from the national inpatient sample (NIS); a US administrative database published under the healthcare utilization project (HCUP) by the Agency for Healthcare Research and Quality ^11^ through a data use agreement.The NIS is the largest publicly available all-payer inpatient care database in the United States, containing data on more than seven million hospital stays from nearly 1,000 U.S. hospitals including non-federal public hospitals and academic medical centers. The data include discharge-level records on demographics, diagnoses, procedures, length of hospital stay (LOS), and resource utilization. The NIS uses the International Classification of Diseases, Clinical Modification 10 (ICD-10) to code for diagnoses and procedures.The data were de-identified, not requiring an IRB approval# The ParticipantsAdult patients with a primary diagnosis of emphysema admitted electively for lung volume reduction (surgical or bronchoscopic). The inclusion and exclusion criteria are aimed at ensuring that the patient was solely admitted for lung volume reduction, and that surgery or endobronchial valve placement was not performed for persistent air leak. Inclusion criteria:* Adults (18 years and older)* Elective admission* Principal diagnosis of emphysema (I10-DX1)* Principal procedure for the admission is lung volume reduction, including bronchoscopic lung volume reduction (I10-PR1) * Principal procedure on the day of admission (procedure on day 0)Exclusion criteria * Non-elective admissions* Presence of persistent air leak diagnosis (ICD10 codes) The sample size consists of 518 adult patients with 253 (49%) patients undergoing bronchoscopic lung volume reduction and 265 patients undergoing surgical lung volume reduction (51%).The study is a retrospective cohort study, two groups of patients with the same baseline disease (emphysema), one group is receiving the intervention of choice (BLVR), and followed during admission for the development of the outcomes **after** the intervention of choice.# The InterventionBronchoscopic lung volume reduction (BLVR) is the intervention of interest. 253 (49%) of the participants out of 518 in my sample have received BLVR and will be compared to surgical lung volume reduction (LVRS), in 265 (51%) of the participants.The intervention (BLVR) will be defined as the first procedure (I10_PR1) received by the participant on day 0 of the elective admission, using ICD-10 procedure codes: 0BH38GZ, 0BH88GZ, 0BH48GZ, 0BH58GZ, 0BH68GZ, 0BH78GZ, 0BH8GZ, 0BH98GZ. The comparison (LVRS) will be defined as the first procedure (I10_PR1) received by the participant on day 0 of the elective admission, using ICD-10 procedure codes: 0BBC0ZZ, 0BBC4ZZ, 0BBD0ZZ, 0BBD4ZZ, 0BBF0ZZ, 0BBF4ZZ, 0BBG0ZZ, 0BBG4ZZ, 0BBH0ZZ, 0BBH4ZZ, 0BBJ0ZZ, 0BBJ4ZZ, 0BBK0ZZ, 0BBK4ZZ, 0BBL0ZZ, 0BBL4ZZ, 0BBM0ZZ, 0BBM4ZZ# The Covariates used in the propensity score The propensity score model used covariates that are related to patient characteristics (age, sex, race, presence of chronic respiratory failure, insurance carrier and classification of the median income for the patient's zip code) and hospital characteristics (hospital size, location, teaching status, hospital ownership and year of the procedure)The covariates were measured and recorded on admission and prior to the procedure.# The Outcomes## Primary outcome: Acute respiratory failure (ARF)Acute respiratory failure is a new (or increased) need for oxygen supplementation in order to maintain adequate oxygenation.It will be defined using ICD-10 diagnosis codes J9600, J9601, J9602, J9620, J9621, J9622, J95821, J95822 and is a binary outcome (Yes/No)## Secondary outcome: Pneumothorax (PTX)Pneumothorax is the abnormal collection of air around the lung leading to lung collapse, and possibly causing respiratory failure or in severe forms, instability and death. This will be defined using ICD-10 diagnosis codes J930, J9311, J9312, J9383, J939, J95811, J95812 and is a binary outcome (Yes/No)# Importing and cleaning the data ## Importing the data The final analytical data was prepared using SAS. This included importing the national inpatient sample (NIS) into a statistical program that can handle the large size and number of original observations (roughly 21 million observations). SAS was used to properly create the variables of interest based on ICD-10 codes. A final dataset with the intended variables was then created for the final analyses. The data was imported in the format of `sas7bdat` using `haven` package. `zap_labels` and `clean_names` were used to remove column labels and standardize names respectively into R format. ```{r, warning=FALSE}blvr <-read_sas("blvr.sas7bdat", NULL) |>zap_label() |>clean_names() ``````{r}glimpse(blvr)```## Cleaning the names and creating factor variables The code below achieved the following steps* Filtering to include adult patients only (Over the age of 18)* Renaming`ebv` to `vol_reduction`, this was intended to be the numeric 0/1 version of the intervention* Mutating character variables as factors and appropriately coding them * Re-leveling certain factors appropriately * The tibble will continue to include `vol_reduction` which is the intervention of choice as a numeric variable (0/1) for further analysis ```{r}blvr <- blvr |>filter(age >17) |>rename(vol_reduction = ebv) |>mutate(sex =fct_recode(factor(female), "Male"="0", "Female"="1"),insurance =fct_recode(factor(pay1), "Medicare"="1", "Medicaid"="2", "Private"="3", "Private"="4", "Private"="6"),insurance =fct_relevel(insurance, "Medicare", "Medicaid", "Private"),race =fct_recode(factor(race), "White"="1", "Black"="2", "Other"="3", "Other"="4", "Other"="5", "Other"="6"),hospital_size =fct_recode(factor(hosp_bedsize), "Small"="1", "Medium"="2", "Large"="3"),teaching =fct_recode(factor(hosp_locteach), "Rural"="1", "Urban nonteaching"="2", "Urban teaching"="3"),region =fct_recode(factor(hosp_region), "Northeast"="1", "Midwest"="2", "South"="3", "West"="4"),ownership =fct_recode(factor(h_contrl), "Government"="1", "Private, non-profit"="2", "Private, for-profit"="3"),income =fct_recode(factor(zipinc_qrtl), "0-25th percentile"="1", "26th to 50th percentile"="2", "51st to 75th percentile"="3", "76th to 100th percentile"="4"),crf =fct_recode(factor(crf), "Yes"="1", "No"="0"),ptx_f =fct_recode(factor(ptx), "Yes"="1", "No"="0"),arf_f =fct_recode(factor(arf), "Yes"="1", "No"="0"),vol_reduction_f =as_factor(volume_reduction),vol_reduction_f =fct_relevel(vol_reduction_f, "LVRS"), year =as_factor(year)) ```## Final analytical tibble The final analytical tibble included the covariates used in building the propensity score, as well as both numeric/ factor versions of the intervention and outcomes ```{r}blvr <- blvr |>select(age, sex, race, crf, income, insurance, hospital_size, teaching, region, ownership, year, vol_reduction, vol_reduction_f, arf, arf_f, ptx, ptx_f)```## Missing variables and Imputation ```{r}blvr |>miss_var_summary()```There were fourteen missing observations in the `race` variable, and nine observations in the `income` variable. These constituted 4.4% of the total data, and assumed to be missing at random (MAR). ```{r}blvr |>tabyl(income)blvr |>tabyl(race)```Given the small sample and the desire to keep as many observations as possible, single imputation using the `mice` package. Fifteen imputed data sets were created and one dataset was later used as the single imputation tibble. ```{r, warning=FALSE}set.seed(123)blvr_mice <-mice(blvr, m =15, printFlag =FALSE)``````{r}blvr <-complete(blvr_mice, 10) |>tibble()```Checking the variables with the missing values after imputation showed preserved ratios in each level in comparison to pre-imputation ```{r}blvr |>tabyl(income)blvr |>tabyl(race)```## Adding an ID row An ID row was added to identify observations if needed ```{r}blvr <- tibble::rowid_to_column(blvr, "ID"); blvr <- blvr |>mutate (ID =as.character(ID))```## Saving final tibbe ```{r}write.csv(blvr, "blvr.csv")```# Table (1) And Codebook ## Final CodebookVariable | Type | Description | Role | -------- | ---- | ----------- | ---- | `age` | Numeric | Age in years | Co-variate | `sex` | Factor (w/2 levels) | Sex of the subject: Male or Female | Covariate | `race` | Factor (w/3 levels) | White, Black, Other | Covariate | `crf` | Factor (w/2 levels) | Presence of chronic respiratory failure: Yes or No | Covariate | `income` | Factor (w/4 levels) | Quartile classification of the estimated median household income of residents in the patient's ZIP Code: 0-25th percentile, 26-50th percentile, 51-75th percentile and 76-00th percentile | Covariate`insurance` | Factor (w/3 levels) | Medicare, Medicaid or Private | Covariate | `hospital_size` | Factor (w/3 levels) | Small, Medium, Large | Covariate |`teaching` | Factor (w/3 levels) | Teaching status of the hospital: Rural, Urban non-teaching and Urban teaching | Covariate |`region` | Factor (w/4 levels) | Region of the hospital: Northeast, Midwest, South and West | Covariate |`ownership` | Factor (w/3 levels) | Ownership of the hospital: Government, Private, non-profit and Private for-profit | Covariate |`year` | Character | Year of the procedure | Covariate | `vol_reduction` | numeric | Type of lung volume reduction: bronchoscopic (1) or surgical (0) | Intervention/Comparator | `vol_reduction_f` | Factor (w/2 levels) | Type of lung volume reduction: bronchoscopic `BLVR` or surgical `LVRS` | Intervention/Comparator | `arf` | Numeric | Presence of acute respiratory failure: 0/1 | Outcome 1 |`arf_f` | Factor (w/2 levels) | Presence of acute respiratory failure: Yes or No | Outcome 1 |`ptx` | Numeric | Presence of pneumothorax: 0/1 | Outcome 1 |`ptx_f` | Factor (w/2 levels) | Presence of pneumothorax: Yes or No | Outcome 2 | ## Creating Table (1) ```{r}blvr.vars <-c("age", "race", "sex", "crf", "income", "insurance","hospital_size", "teaching", "region", "ownership", "year")blvr.reduction <-c("vol_reduction_f")blvr_table1 <-CreateTableOne(data = blvr, vars = blvr.vars,strata = blvr.reduction,test =FALSE)``````{r}print(blvr_table1) |>kable() |>kable_material()```Based on table 1, patients receiving bronchoscopic lung volume reduction (BLVR) were noted to be older (mean age of 67 years), with a higher percentage of females (43.5%), likely due to a preference of a less invasive approach compared to the surgical group. BLVR patients also tended to have a higher percentage of chronic respiratory failure, likely due to the ease of valve removal should the patient's baseline respiratory failure worsen after procedure. The patient population was more white (90%), with medicare as the primary insurance provider (77.5%) which is likely due to the age bracket in this group (65 years and older). As expected both groups had a similar distribution of median income per zip code, and majority of these procedures were done in large, non for profit and urban teaching hospitals. Lung volume reduction surgery number distribution was steady over the years included (roughly around 30%), however bronchoscopic lung volume reduction numbers significantly increased over the years (starting in 2019 and 2020)# Unadjusted Outcomes ## Acute Respiratory Failure `arf` ```{r}ggplot(blvr, aes(x = vol_reduction_f, fill = arf_f)) +geom_bar(aes(y = (after_stat(count))/sum(after_stat(count))), position ="dodge") +scale_y_continuous(labels = scales::percent) +scale_fill_viridis_d(option ="turbo", alpha =0.5) +labs(title ="Respiratory Failure", x ="Volume Reduction", y ="% of Patients", fill ="Resp. Failure")``````{r}blvr |>tabyl(arf_f, vol_reduction_f) |>kable() |>kable_material()```The visualization above suggested that BLVR patients tended to suffer a higher percentage of acute respiratory failureA logistic regression model was fit to assess the odds of respiratory failure based on the type of volume reduction received, then the coefficient was exponentiated to obtain an odds ratio ```{r}unadjusted_arf <-glm(arf ~ vol_reduction, data=blvr, family=binomial)``````{r}kable(tidy(unadjusted_arf, conf.int =TRUE, conf.level =0.95, exponentiate =TRUE), digits =2) |>kable_paper()```If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 1.65 odds of acute respiratory failure (95% CI (1,2.74)). This confidence interval contained 1, thus it was difficult to reliably say that BLVR would meaningfully lead to a higher odds of acute respiratory failure Overall, it seems that the patients receiving BLVR tended to have higher odds of respiratory failure, even though the odds ratio's interval contained 1, this odds ratio shows an interesting trend ## Pneumothorax `ptx` ```{r}ggplot(blvr, aes(x = vol_reduction_f, fill = ptx_f)) +geom_bar(aes(y = (after_stat(count))/sum(after_stat(count))), position ="dodge") +scale_y_continuous(labels = scales::percent) +scale_fill_viridis_d(option ="turbo", alpha =0.5) +labs(title ="Pneumothorax", x ="Volume Reduction", y ="% of Patients", fill ="PTX")``````{r}blvr |>tabyl(ptx_f, vol_reduction_f) |>kable() |>kable_material()```The visualization above suggested that BLVR patients tended to suffer a lower percentage of pneumothoraxA logistic regression model was fit to assess the odds of pneumothorax based on the type of volume reduction received, then the coefficient was exponentiated to obtain an odds ratio ```{r}unadjusted_ptx <-glm(ptx ~ vol_reduction, data=blvr, family=binomial)``````{r}kable(tidy(unadjusted_ptx, conf.int =TRUE, conf.level =0.95, exponentiate =TRUE), digits =2) |>kable_paper()```If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR had a 0.4 odds of pneumothorax (95% CI (0.27,0.58)), which is lower odds than the subject receiving LVRS, which is plausible as LVRS is more invasive. # Propensity Score ## Fitting the model A logistic regression model was used to estimate the propensity score. Having BLVR/LVRS was modeled as the outcome, with patient and hospital characteristics as independent variables. Patient characteristics used in the propensity score included `age`, `sex`, `race`, `crf`, `income` and `insurance`. Hospital characteristics included `hospital_size`, `teaching` status, `region`, `ownership` and `year` of the procedure. The raw and linear propensity score values were saved to the `blvr` tibble. ```{r}psmodel <-glm(vol_reduction ~ age + sex + race + crf + income + insurance + hospital_size + teaching + region + ownership + year,family=binomial, data=blvr)summary(psmodel)``````{r}blvr$ps <- psmodel$fittedblvr$linps <- psmodel$linear.predictors```## Describing the overlap ```{r}ggplot(blvr, aes(x = vol_reduction_f, y = ps, color = vol_reduction_f)) +geom_boxplot() +geom_jitter(width =0.1, alpha =0.3) +scale_color_viridis_d(option ="plasma", end =0.5) +labs(title ="Boxplot & Jitter", subtitle ="Distribution Of the Raw PS", x ="Type of Procedure", y ="Propensity Scores",caption ="* Line indicates sample median", col ="Procedure")``````{r}ggplot(blvr, aes(x = linps, fill = vol_reduction_f)) +geom_density(alpha =0.3) +scale_fill_viridis_d(option ="turbo") +labs(title ="Density Plot", subtitle ="For the Linear Propensity Score", x ="Linps", y ="Density",fill ="Volume Reduction")``````{r}blvr |>group_by(vol_reduction_f) |>skim_without_charts(ps, linps)```The overlap between the propensity scores in intervention groups was assessed numerically and graphically. The boxplot and the density plots show some overlap in the propensity scores between the intervention and the treatment.The LVRS group had a lower propensity score mean, with problematic values falling below 0.01. The BLVR group also had some problematic high values over 0.99. Troubleshooting will be addressed at a later step ## Unadjusted Rubin's Rules ```{r}unadjusted_rubin1 <-abs(100*(mean(blvr$linps[blvr$vol_reduction==1]) -mean(blvr$linps[blvr$vol_reduction==0])) /sd(blvr$linps))unadjusted_rubin2 <-with(blvr, var(linps[vol_reduction ==1]) /var(linps[vol_reduction ==0]))``````{r}unadjusted_rubin1; unadjusted_rubin2```Rubin's rule 1 value was higher than 50%, this meant that running an unadjusted model could not be justified without applying the propensity scores to account for selection bias. Rubin's rule 2 was on the lower bound of acceptable value as well, as the variance ratio of the linear propensity score in the BLVR (intervention) group to the variance of the linear propensity score in the LVRS (control) group should be close to 1 (and ideally between 0.8 and 1.25)## Troubleshooting the Extreme Propensity Score Values ```{r}blvr |>count(ps >=0.99); blvr |>count(ps <0.01)```There were twenty eight observations with propensity score values below 0.01 and four observations with propensity score values above 0.99. The extreme lows seemed to be present in the LVRS group while the extreme highs were present in the BLVR group ```{r}mosaic::favstats(ps ~ vol_reduction_f, data = blvr)```The following strategies were attempted in order to resolve the above issue:- Colinearity check using `vif` didn't reveal problematic colinearity between the variables used in the propensity score model ```{r}car::vif(psmodel)```- Re-constructing the propensity score model, adding one covariate at a time, while checking which variable drives the propensity score values below 0.1 or above 0.99 showed that each of the following covariates `race`, `year`, `sex`, `ownership` and `crf` could independently drop the propensity score below 0.01 or drive it up above 0.99. Excluding these covariates form the model would be problematic in itself - Collapsing multicategorical variables was attempted, these variables included `year`, `race`, `ownership`, `insurance`, `hospital_size`, `teaching` and `region`. It did not result in the improvement of the propensity score even when each covariate was added individually after collapsing. Filtering obaservations with problematic propensity scores proved to be the most reasonable approach in this situation. This resulted in a sample with 237 observations on the control intervention (LVRS) and 249 observation on the studied intervention (BLVR). The resulting dataset `blvr_matching` was used for propensity score matching analyses```{r}blvr_matching <- blvr |>filter(ps <=0.990& ps >=0.010)``````{r}mosaic::favstats(ps ~ vol_reduction_f, data = blvr_matching)```A density plot was created to describe the overlap of the linear propensity scores between the two intervention groups after balancing the dataset, the overlap is described graphically in the density plot below.```{r}ggplot(blvr_matching, aes(x = linps, fill = vol_reduction_f)) +geom_density(alpha =0.3) +scale_fill_viridis_d(option ="turbo") +labs(title ="Density Plot", subtitle ="For the Linear Propensity Score", x ="Linps", y ="Density",fill ="Volume Reduction")```# Approach (1): Matching {#sec-matching}Propensity score matching was conducted using a 1:1 caliper matching on the linear propensity score without replacement from the `Matching` package. A default caliper value of 0.2 was selected. The type of matching was chosen based on the similar sample size between the intervention and control groups. ```{r}X <- blvr_matching$linpsTr <-as.logical(blvr_matching$vol_reduction)match <-Match(Tr = Tr, X = X, M =1, caliper =0.2,estimand ="ATT", replace =FALSE, ties =FALSE)summary(match)```One hundred and seven BLVR patients were matched to one hundred and seven LVRS patients. The degree of covariate imbalance was assessed before and after propensity score matching for each covariate as well as the linear and raw propensity scores ```{r}covs <- blvr_matching |>select(age, sex, race, crf, insurance, income, hospital_size, region, teaching, ownership, year, ps, linps)b <-bal.tab(match,treat = blvr_matching$vol_reduction,covs = covs,stats =c("m", "v"),quick =FALSE, un =TRUE, disp.v.ratio =TRUE)b``````{r, message=FALSE}bal.plot(match,treat = blvr_matching$vol_reduction,covs = covs,var.name ="linps",which ="both",sample.names =c("Unmatched Sample", "Matched Sample")) +scale_fill_viridis_d(option ="turbo") +labs(title ="Distributional Balance for linear PS", subtitle ="1:1 Caliper Matching", x ="Linear PS") ```A Love plot of standardized differences before and after matching was obtained to evaluate matching success, showing improvement over the unmatched sample ```{r}love.plot(b, threshold = .1, size =3,var.order ="unadjusted",colors =c("hotpink3", "purple4"),stars ="raw") +labs(title ="Love Plot of Standardized Differences", subtitle ="1:1 Caliper Matching")```Rubin's first and second rules were assessed after matching, showing improvement over the unadjusted Rubin's rules.```{r}matches <-factor(rep(match$index.treated, 2))blvr.matchedsample <-cbind(matches, blvr_matching[c(match$index.control, match$index.treated),])``````{r}matched_rubin1 <-with(blvr.matchedsample,abs(100*(mean(linps[vol_reduction==1])-mean(linps[vol_reduction==0]))/sd(linps)))matched_rubin2 <-with(blvr.matchedsample, var(linps[vol_reduction==1])/var(linps[vol_reduction==0]))``````{r}matched_rubin1; matched_rubin2```Rubin's rule 1 was less than 10%, and Rubin's rule 2 was around 1 which was ideal. # Outcomes after Matching ## Acute respiratory failure `arf````{r}matched_arf <-clogit(arf ~ vol_reduction +strata(matches), data=blvr.matchedsample)``````{r}kable(tidy(matched_arf, conf.int =TRUE, conf.level =0.95, exponentiate =TRUE), digits =2) |>kable_material()```If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 1.73 odds of acute respiratory failure (95% CI (0.82,3.63)). The subject receiving BLVR tended to have a higher odds ratio of respiratory failue in comparison to patients receiving LVRS, however this confidence interval crosses 1## Pneumothorax `ptx````{r}matched_ptx <-clogit(ptx ~ vol_reduction +strata(matches), data=blvr.matchedsample)``````{r}kable(tidy(matched_ptx, exponentiate =TRUE, conf.level =0.95, conf.int =TRUE), digits =2) |>kable_material()```If we have two subjects, one who received LVRS and one who received BLVR, then the subject receiving BLVR has a 0.68 odds of pneumothorax (95% CI (0.38,1.22)). Patients receiving BLVR seemed to have a higher odds ratio of pneumothorax in comparison to patients receiving BLVR, however this confidence interval contained 1, the results were not statistically detectable and a stability analysis instead of a sensitivity analysis was carried out. ### Stability Analysis- Trying different imputed datasets (generated through `mice` package) resulted in similar outcome estimates (odds ratio)- The matching strategies described below were attempted. None of these strategies provided an improvement in the effective sample size, propensity score distributional balance, standardized differences or Rubin's rules in comparison to the matching strategy used in this analysis (1:1 matching with a caliper)Match | Trtmnt | Ctrl | Distributional Balance | Std Differences | Rubin's 1 | Rubin's 2 | ----- | ------ | ---- | -----------------------| --------------- | --------- | --------- | `1:1 Greedy W replacement` | 249 | 75 | Better overlap, not ideal | better but dispersed | 8.8 | 1.27 | `Genetic matching` | 273 | 78 | Better overlap, not ideal | Dispersed | 8.8 | 1.26 | `1:1 Neartest neighbor matching W Replacement` | 249 | 79 | Better overlap, not ideal | better but dispersed (> 0.1) | 8 | 1.21 |- The following steps attempted comparing the odds ratios obtained through matching to the odds ratios after ATT weighting and a double robust approach (by adjusting for the linear propensity score in the regression model) # Approach (2): Propensity Score Weighting {#sec-weighting}Given the initial difference in sample sizes across the two volume reduction cohorts, propensity score matching has left a selection of the patients unmatched. Weighting by the inverse propensity score (average of the treatment effect on the treated) was next carried out using the initial database `blvr` (the database that still contained propensity scores below 0.01 and above 0.99)```{r}blvr$wts1 <-ifelse(blvr$vol_reduction==1, 1, blvr$ps/(1-blvr$ps))```Then the resulting ATT (average treatment effect on the treated) weights were plotted in the figure below, where the treatment group got the same weight of 1, and the control group had varying weights ```{r}ggplot(blvr, aes(x = ps, y = wts1, color = vol_reduction_f)) +geom_point(size =2, alpha =0.4) +facet_wrap(~ vol_reduction_f) +labs(x ="Estimated Propensity for Volume Reduction",y ="ATT weights for blvr",title ="ATT weighting structure",subtitle ="In the BLVR dataset",col ="Volume Reduction")```The balance after weighting was assessed using `dx.wts` from the `twang` package. Prior to weighting, the control sample was 237 participants, and after weighting, the effective sample size was 54 participants, representing a sample size that would be required to achieve the same level of precision if that sample were a simple random sample```{r, warning=FALSE}blvr_df <-data.frame(blvr)covlist <-c("age", "sex", "race", "crf", "insurance", "income", "hospital_size", "teaching", "region", "ownership", "year", "ps", "linps") bal.wts1 <-dx.wts(x=blvr_df$wts1, data=blvr_df, vars=covlist, treat.var="vol_reduction", estimand="ATT")bal.wts1```The standardized differences were compared before and after weighting using a `Love.plot`. The weighting process showed a more satisfactory balance of standardized differences in comparison to the unweighted standardized differences. ```{r}bal.table(bal.wts1)``````{r}bal.before.wts1 <-bal.table(bal.wts1)[1]bal.after.wts1 <-bal.table(bal.wts1)[2]balance.att.weights <-data.frame(names =rownames(bal.before.wts1$unw), pre.weighting =100*bal.before.wts1$unw$std.eff.sz, ATT.weighted =100*bal.after.wts1[[1]]$std.eff.sz)balance.att.weights <-gather(balance.att.weights, timing, szd, 2:3)balance.att.weights``````{r}ggplot(balance.att.weights, aes(x = szd, y =reorder(names, szd), color = timing)) +geom_point(size =2, alpha =0.5) +geom_vline(xintercept =0) +geom_vline(xintercept =c(-10,10), linetype ="dashed", col ="blue") +labs(x ="Standardized Differences", y ="", title ="Love Plot of Standarized Diffs in The BLVR Data",subtitle ="Before and After ATT Weighting") ```Rubin's rules were obtained after weighting. Rubin's rule 1 was obtained by using the standardized effect size after weighting for the linear propensity score (0.28). When multiplied by 100, the result was 28% indicating that Rubin's rule 1 was satisfied (although not perfectly). Rubin's rule 2 was calculated by dividing the square root of treatment group standard deviation by the square root of control group standard deviation, resulting in a value of 1.54. Rubin's rule 2 was also satisfied, although less perfectly. # Outcomes after Weighting ## Acute Respiratory Failure `arf````{r}blvr.design <-svydesign(ids=~1, weights=~wts1, data=blvr)``````{r}weighted_arf <-svyglm(arf ~ vol_reduction, design=blvr.design, family=quasibinomial())``````{r}kable(tidy(weighted_arf, conf.int =TRUE, conf.level =0.95, exponentiate =TRUE), digits =2) |>kable_material() ```The estimate of treatment effect on acute respiratory failure was obtained by a regression model using survey weighting and a quasi binomial design.The weighted odds ratio estimate for the BLVR group was 1.68, with a 95% confidence interval (0.83,3.4). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having acute respiratory failure for the person receiving BLVR was higher and 1.68 the odds of respiratory failure compared to the person receiving LVRS. It should be noted that this confidence interveal crosses 1 ## Pneumothorax `ptx````{r}weighted_ptx <-svyglm(ptx ~ vol_reduction, design=blvr.design, family=quasibinomial())``````{r}kable(tidy(weighted_ptx, conf.int =TRUE, conf.level =0.95, exponentiate =TRUE), digits =2) |>kable_material() ```The weighted odds ratio estimate for the `BLVR group was 0.36, with a 95% confidence interval (0.19,0.67). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having pneumothorax for the person receiving BLVR was lower and was 0.36 the odds of pneumothorax of the person receiving LVRS. The confidence interval did not cross 1# Apprach (3): A Double-Robust Estimate A logistic regression for the outcomes of acute respiratory failure and pneumothorax was fitted with a regression adjustment for the linear propensity score to obtain a doubly-robust estimate and account for residual imbalance after weighting ```{r}robust.design <-svydesign(ids=~1, weights=~wts1, data=blvr)``````{r}robust_arf <-svyglm(arf ~ vol_reduction + linps, design=robust.design, family=quasibinomial())``````{R}kable(tidy(robust_arf, exponentiate =TRUE, conf.int =TRUE, conf.level =0.95), digits =2) |>kable_material()```The doubly-robust odds ratio estimate for the BLVR group was 1.89, with a 95% confidence interval (0.93,3.87). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having acute respiratory failure for the person receiving BLVR was higher and 1.89 the odds of acute respiratory failure compared to the person receiving LVRS. The confidence interval here crosses 1 ## Pneumothorax `ptx````{r}robust_ptx <-svyglm(ptx ~ vol_reduction + linps, design=robust.design, family=quasibinomial())```And then the Odds Ratio with the 95% confidence interval ```{R}kable(tidy(robust_ptx, conf.int =TRUE, exponentiate =TRUE, conf.level =0.95), digits =2) |>kable_material()```The doubly-robust odds ratio estimate for the `BLVR group was 0.38, with a 95% confidence interval (0.2,0.75). If we had two subjects, one who had BLVR and one who had LVRS, then the odds of having pneumothorax for the person receiving BLVR was lower and was 0.38 the odds of pneumothorax of the person receiving LVRS. The confidence interval did not cross 1# Final comparison of Outcome Estimates ## Acute Respiratory Failure `arf` {#sec-arf}```{r}arf_blvr <-tibble(analysis =c("Unadjusted", "Matched", "ATT Weighted", "Double-Robust"),estimate =c(1.65, 1.73, 1.68, 1.89),conf.low =c(1, 0.82, 0.83, 0.93),conf.high =c(2.74, 3.63, 3.4, 3.87))ggplot(arf_blvr, aes(x = analysis, y = estimate)) +geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width =0.5, color ="orange3") +geom_label(aes(label = estimate), size =6, color ="white", fill ="lightblue3") +labs(title ="Comparison Of Odds Ratio Estimates", subtitle ="Between Different Adjustment Approaches In Respiratory Failure", x ="Method of PS Adjustment", y ="Odds Ratio")```Method | Estimate | Lower 95% CI | upper 95% CI |------ | -------- | ------------ | ------------ | Unadjusted | 1.65 | 1.00 | 2.74 |Matched | 1.73 | 0.82 | 3.63 | Weighted | 1.68 | 0.83 | 3.4 | Double-Robust | 1.89 | 0.93 | 3.87 |Based on the comparison above, the odds ratio of acute respiratory failure didn't substantially change from its estimate prior to propensity score adjustment. All of the confidence intervals contain 1.## Pneumothorax `ptx` {#sec-ptx}```{r}ptx_blvr <-tibble(analysis =c("Unadjusted", "Matched", "ATT Weighted", "Double-Robust"),estimate =c(0.4, 0.68, 0.36, 0.38),conf.low =c(0.27, 0.38, 0.19, 0.2),conf.high =c(0.58, 1.22, 0.67, 0.75))ggplot(ptx_blvr, aes(x = analysis, y = estimate)) +geom_errorbar(aes(ymax = conf.high, ymin = conf.low), width =0.5, color ="purple3") +geom_label(aes(label = estimate), size =6, color ="white", fill ="lightpink3") +labs(title ="Comparison Of Odds Ratio Estimates", subtitle ="Between Different Adjustment Approaches In Pneumothorax", x ="Method of PS Adjustment", y ="Odds Ratio")```Method | Estimate | Lower 95% CI | upper 95% CI |------ | -------- | ------------ | ------------ | Unadjusted | 0.4 | 0.27 | 0.58 | Matched | 0.68 | 0.38 | 1.22 | Weighted | 0.36 | 0.19 | 0.67 | Double-Robust | 0.38 | 0.2 | 0.75 | The comparison above showed that the odds ratio for pneumothorax also didn't substantially change after propensity score adjustment, reflecting a plausible conclusion that patients undergoing lung volume reduction surgery tended to have higher odds of pneumothorax based on the nature of invasive surgical intervention on the lung. # Consclusions ## Answering The Research Questions > Is bronchoscopic lung volume reduction associated with higher odds of developing post-procedural acute respiratory failure in adult patients with emphysema compared to lung volume reduction surgery?Yes. Based on the odds ratio obtained prior to propensity score adjustment and after, the patients undergoing bronchoscopic lung volume reduction surgery tended to have higher odds of developing acute respiratory failure, granted that the confidence interval for this estimate always crossed 1 It was noted that despite propensity score matching, weighting and a double robust adjustment (@sec-arf), the odds ratio of acute respiratory failure did not show a meaningful change from prior to adjustment. The small fluctuations wouldn't have significant clinical implication. Even though the confidence intervals crossed 1, this analysis showed an interesting trend in the odds of respiratory failure. A possible explanation could be that there was a higher percentage of patients with chronic respiratory failure undergoing bronchoscopic lung volume reduction in comparison to patients undergoing surgical lung volume reduction, so the patient's lungs have already been compromised and further interventions > Is bronchoscopic lung volume reduction associated with a higher odds of developing post-procedural pneumothorax in adult patients with emphysema compared to lung volume reduction surgery?No. BLVR was associated with lower odds of pneumothorax (@sec-ptx) despite propensity score adjustment methods, which was rather plausible, as the surgical procedure was more invasive and expected to cause a higher incidence of pneumothorax. The presence of pneumothorax is well described in the literature, and could indicate a favorable response after lung volume reduction ## Future Directions```{r}ggplot(blvr, aes(x = year, fill = vol_reduction_f)) +geom_bar(position ="dodge", col ="white") +scale_fill_viridis_d(option ="turbo", alpha =0.5) +labs(title ="Lung Volume Reduction", subtitle ="Distribution By Year", x ="Volume Reduction", y ="# of Patients", fill ="Volume Reduction") ```Bronchoscopic lung volume reduction with endobronchial valves has largely replaced lung volume reduction surgery as a possible solution to improve excessive dyspnea. Patient selection is a key element in developing post procedural complications that could be life threatening. Understanding the potential complications that may arise in real-world settings outside controlled clinical trials becomes increasingly crucial as the number of procedures grows and their availability expands beyond large referral centers. This knowledge aids physicians in better patient selection for the procedure, enabling them to implement effective action plans. For instance, the introduction of a "Pneumothorax Travel Cart" containing procedural kits for pneumothorax management has significantly improved time to management of pneumothorax after bronchoscopic lung volume reduction. By trying to get a sense of the complication rates nationally, this information could be used to obtain further data from endobronchial valve device companies. These data could include more patient covariates that play a role in procedural selection, and provide a better characterization of the complications encountered. ## Learning Curve By doing this project, I got a chance to apply all the matching methods (@sec-matching) we learned over the span of this course. While I was slightly distraught at the beginning when I ended up with a sample with more treated patients than controls, I used it as an opportunity to go over all the plausible methods to properly troubleshoot and be prepared for similar scenarios in the future beyond the security of this class. # Acknowledgments - Dr. Marc Georgi, a gastroenterologist at Inova Fairfax Hospital, for simplifying the use of the Healthcare Utilization Project databases (the national inpatient sample and the national readmission database). - Dr. Thomas Love, for making statistics incredibly enjoyable!# References 1. Global strategy for the diagnosis, management and prevention of chronic obstructive pulmonary disease. Global Initiative for Chronic Obstructive Lung Disease; 20212. Fishman A, Martinez F, Naunheim K, et al. A randomized trial comparing lung-volume-reduction surgery with medical therapy for severe emphysema. N Engl J Med. 2003;348(21):2059-20733. Fernandez-Bussy S, Labarca G, Herth FJF. Bronchoscopic Lung Volume Reduction in Patients with Severe Emphysema. Semin Respir Crit Care Med. 2018;39(6):685-692. 4. Criner GJ, Sue R, Wright S, et al. A Multicenter Randomized Controlled Trial of Zephyr Endobronchial Valve Treatment in Heterogeneous Emphysema (LIBERATE). Am J Respir Crit Care Med. 2018;198(9):1151-1164.5. Criner GJ, Delage A, Voelker K, et al. Improving Lung Function in Severe Heterogenous Emphysema with the Spiration Valve System (EMPROVE). A Multicenter, Open-Label Randomized Controlled Clinical Trial. Am J Respir Crit Care Med. 2019;200(11):1354-13626. Zantah M, Gangemi AJ, Criner GJ. Bronchoscopic lung volume reduction: status quo. Annals of translational medicine 2020; 8:1469.7. Davey C, Zoumot Z, Jordan S, McNulty WH, Carr DH, Hind MD, Hansell DM, Rubens MB, Banya W, Polkey MI, et al. Bronchoscopic lung volume reduction with endobronchial valves for patients with heterogeneous emphysema and intact interlobar fissures (the BeLieVeR-HIFi study): a randomised controlled trial. Lancet 2015;386:1066–10738. Valipour A, Slebos DJ, Herth F, Darwiche K, Wagner M, Ficker JH, Petermann C, Hubner RH, Stanzel F, Eberhardt R; IMPACT Study Team. Endobronchial valve therapy in patients with homogeneous emphysema: results from the IMPACT study. Am J Respir Crit Care Med 2016;194:1073–1082.9. Fiorelli A, D'Andrilli A, Bezzi M, et al. Complications related to endoscopic lung volume reduction for emphysema with endobronchial valves: results of a multicenter study. J Thorac Dis. 2018;10(Suppl 27):S3315-S3325. doi:10.21037/jtd.2018.06.6910. Posthuma R, Vaes AW, Walraven KHM, et al. Implementation of Bronchoscopic Lung Volume Reduction Using One-Way Endobronchial Valves: A Retrospective Single-Centre Cohort Study. Respiration. 2022;101(5):476-484. doi:10.1159/00052088511. Healthcare Cost and Utilization Project. HCUP National Inpatient Sample (NIS). Agency for Healthcare Research and Quality; 2019-2020.# Session Information```{r}xfun::session_info()```